function y = lik(par,Ndem,Nrep,N)
global p_values gammapar
theta=par(1); alpha=par(2);
disp(['Coefficients: theta=',num2str(theta),'; alpha=', num2str(alpha)])
options = optimset('TolFun',1e-6,'TolX',1e-6);

y=zeros(size(N));
p_values=zeros(size(N));
for i=1:1:size(N,1)
    g=@(x,y) exp(new_gammaln(N(i)+1)-new_gammaln(x+1)-new_gammaln(N(i)-x+1)+betaln(x+theta,N(i)-x+theta)-betaln(theta,theta)+new_gammaln(x+1)-new_gammaln(Ndem(i)+1)-new_gammaln(x-Ndem(i)+1)+Ndem(i).*log(max(min(y.*(x./(N(i)-x)).^gammapar,1),0))+(x-Ndem(i)).*log(1-max(min(y.*(x./(N(i)-x)).^gammapar,1),0))+new_gammaln((N(i)-x)+1)-new_gammaln(Nrep(i)+1)-new_gammaln((N(i)-x)-Nrep(i)+1)+Nrep(i).*log(max(min(alpha.*y.*((N(i)-x)./x).^gammapar,1),0))+((N(i)-x)-Nrep(i)).*log(1-max(min(alpha.*y.*((N(i)-x)./x).^gammapar,1),0)));
    t=@(y) -(quad(@(x)g(x,y),Ndem(i),N(i)-Nrep(i),1e-12));    
    p=1;
    step=2*Ndem(i)/N(i);
    while p>step-0.0001 & step<0.15
    step=step+0.001;
    p=fminbnd(t,1e-20,step,options);
    end
    if step>.14
        error
    end
    if mod(i,100)==0
    fprintf('%d%% done. ', floor((i/size(N,1))*100));
    end
    p_values(i)=p;
    y(i)=quad(@(x)g(x,p),Ndem(i),N(i)-Nrep(i),1e-12);
end
nonzero=find(y>0); zero=mean(y==0)
y=-sum(log(y(nonzero)));
end